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Abstract 

The diffusion process near low order synchro-betatron resonances driven by beam- 
beam interactions at a crossing angle is investigated. Macroscopic observables such as 
beam emittance, lifetime and beam profiles are calculated. These are followed with de- 
tailed studies of microscopic quantities such as the evolution of the variance at several 
initial transverse amplitudes and single particle probability distribution functions. We 
present evidence to show that the observed diffusion is anomalous and the dynamics 
follows a non-Markovian continuous time random walk process. We derive a modified 
master equation to replace the Chapman-Kolmogorov equation in action-angle space 
and a fractional diffusion equation to describe the density evolution for this class of 
processes. 

1 Introduction 

Diffusion of particle beams due to nonlinear fields is often a major source of emittance 
growth and beam loss in an accelerator. Measurements of diffusion coefficients have been 
reported from several hadron accelerators [HJ E E]. The diffusion equation was also used to 
explain the change in beam lifetime following the failure of a separator during a Tevatron 
store [4j. In collision mode the beam-beam interactions are usually the dominant nonlin- 
earity. Diffusion coefficients in the absence of low order resonances have been calculated 
for head-on interactions [5] and for long-range interactions O. Diffusion due to nonlinear 
resonances is more complex and the study of this phenomenon has a long history, see e.g 
[HI [9l PTOl tTTJ . Resonances when modulated, either by dynamical effects such as synchro- 
betatron coupling or due to ripple in magnet currents, can sweep across phase space and 
transport particles to large amplitudes ||12|. 

In this article we will study the nature of the diffusion process due to synchro-betatron 
resonances driven by beam-beam interactions with a crossing angle. This was first inves- 
tigated at the DORIS collider [fT3l and has since been observed at other colliders. Our 



1 



aim is to establish the correct statistical mechanical model that describes the evolution of 
the beam density. We examine the possibility that the diffusion process is anomalous with 
detailed tracking simulations and derive a master equation and a related fractional diffu- 
sion equation that may describe the transport process. A preliminary version of this study 
was reported in [|T4l . An example of anomalous diffusion observed in particle beams as a 
consequence of rf phase modulation was reported in [15]. Anomalous diffusion processes 
have been reported in several areas of physics including plasma turbulence lfT6l . and in the 
motion of laser cooled atoms on a lattice [fTTl . 

2 Synchro-betatron resonances due to crossing angles 

Synchro-betatron resonances (SBRs) due to beam-beam interactions at a crossing angle 
are convenient to study resonantly driven amplitude growth for several reasons. At large 
amplitudes, the non-linear force vanishes, hence particle excursions do not go to arbitrarily 
large amplitudes which is not the case for resonances due to multipole nonlinearities. This 
removes numerical instabilities and also allows the entire beam to be probed for the particle 
dynamics. Another advantage is that the resonances can be studied in one transverse plane 
since these resonances are driven by energy pumped from the longitudinal plane to the 
transverse plane with very little impact on the longitudinal dynamics. 

When beams collide at an angle, the transverse distance of a test particle from the center 
of the opposing bunch depends on the longitudinal position of the particle. Consequently 
synchrotron oscillations of the particle couple to the transverse beam-beam force leading 
to excitation of synchro-betatron resonances. Since the beam-beam force goes to zero at 
large transverse separations, the effects of these resonances are experienced by particles 
only within a certain range of transverse amplitudes. 

For simplicity, we choose the resonances to be in only one transverse plane, here the 
horizontal plane. In order to observe effects over relatively short computation times, we 
choose low order resonances. The tunes we choose are unrealistic for operating colliders 
but it is likely that the dynamics near high order resonances is similar but occurs over a 
longer time scale. 

Linear motion and the beam-beam interactions can be described by the equations of 
motion resulting from the Hamiltonian 

N IP 

H = v x J x + v y J y + v s J s + £ Ui (x, y,s)8p(Q- $) (1) 

j 

where (v x ,V y ,V s ) are the tunes, and (J x ,Jy,Js) are the actions. U(x,y,s) is the beam-beam 
potential, dp is the periodic delta function, is the azimuthal coordinate and the sum 
extends over the number Njp of interaction points. Assuming Gaussian distributions in all 
three planes, crossing angles of (20*, 20y) in the horizontal and vertical planes respectively, 



2 



the beam-beam potential for colliding proton bunches can be written as 



rrl \ Nbr P t 
U{x,y,s) = — — / 

Yd Jo 



N b r p [°° dq 



Yp Jo [(2ffi + q)[2o* + q)]W 

(x + ssinlfa) 2 (y + ssm2(j) y ) : 

where Nf, is the bunch intensity of the opposing bunch, r p is the classical proton radius, 
Y p is the proton energy in units of its rest mass and o x , o y are the rms beams sizes of the 
opposing beam at the interaction point (IP). The potential can be expanded as a Fourier 
series 

U(x,y,s)= U m x ,m y ,m s ^P[K m x¥x + fny¥y + m s¥s-P^)] (3) 

m x ,m y ,m s ,p 

This potential can excite synchro-betatron resonances given by the resonance condition 
m x V x + niyVy + m s v s = p where (m Xl m yi m Sl p) are integers. It can be shown from the 
structure of the Fourier harmonics U m ^ m ms that they are non-zero only when the sum 
m x + m y + m s is even. The Fourier harmonics can also be used to calculate the tune shifts 
with amplitude and the resonance driving terms, as was done in lfT8l . As one example, we 
write down the zero transverse amplitude tune shift for round beams. This tune shift now 
depends on the longitudinal oscillation amplitude a s o s as 

Av x (a x = 0,a y = 0,a s ) = ^[I (t) - /i(t)(1 + \^-)\ (4) 

and a similar expression for Av y . Here £ = ~N r p l is the usual beam-beam param- 
eter, (a x <J x ,a y Oy) are the transverse amplitudes of the particle, IqJ\ are modified Bessel 
functions and the other dimensionless parameters are 

z = -aAh x +ht), h x = — sin2^ x , h y = — sin20 v 

4 • y G X Oy 

As a consequence, only those zero transverse amplitude particles with zero longitudinal 
amplitude a s experience the full beam-beam tune shift £ . Particles with non-zero amplitude 
a s experience a smaller tune shift. 

Since the LHC employs crossing angles in its collision scheme, we will use the LHC 
beam parameters in the simulations reported here. As in the LHC, the crossing angle is in 
the horizontal plane at one IP and in the vertical plane at the second IP. We consider reso- 
nances excited in the horizontal plane only, so they are of the form m x v x + m s v s = p with 
m x +m s even. In our model the only sources of tune spread are the beam-beam interactions. 
These interactions between protons lowers the betatron tunes at small amplitudes. We 
choose the large amplitude tunes, i.e. the tunes with only the linear lattice, to satisfy one of 
the SBR resonance conditions. Having chosen a particular resonance m x V x + mu s V s = p to 
be satisfied by the bare lattice tunes, the tunes inside the bunch are determined by the beam- 
beam parameter the synchrotron tune V s and the amplitudes (a x ,a y ,a s ) of the particle. 
The nominal LHC horizontal tune is 0.31 at collision, so we searched among the following 
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Table 1: Table of basic parameters in simulation model. Resonance I is 2(3v x — 2v s ) = 2, 
resonance II is 4v x — 2v s = 1 . 



Beam parameter 


Value 


bnergy UeVJ 


7.0 


Bunch Intensity 


l.lxlO 11 


Ox, o y [jum] 


16.6, 16.6 


CTj [cm] 


7.5 


Rf voltage [MV] 


16 


Crossing angles [/irad] 


300 


Beam-beam parameter 


0.0034 


Resonance I: (v x ,Vy) 


0.3353, 0.32 


Resonance II: (v x , Vy) 


0.2514, 0.32 



resonances: 3v x ± V s = 1, 2(3v x ±2v s ) — 2 as well as 2(4v x ± V s ) = 2 and 4v x ±2v s = 1 to 
find those that cause large growth of the emittance and beam tails. Given that the betatron 
tune spread from head-on beam-beam interactions is about 0.007 and the small amplitude 
synchrotron tune is ~0.002, the choices 2(3v x — 2v s ) = 2 and 4v x — 2v s = 1 had the greatest 
impact on the beam. With these choices, low amplitude particles are resonant with the third 
and fourth order betatron resonances respectively, and the synchrotron oscillations modu- 
late these resonances leading to large amplitude growth. The other resonances are resonant 
at larger amplitudes and consequently have a smaller impact on the bunch. The bare lattice 
(which become the large amplitude) betatron tunes corresponding to these resonances are 
shown in Table [Q Some of these parameters may be slightly different from the present 
LHC design values, e.g the LHC design value of the crossing angle is 285/irad. 

3 Simulations of beam variables 

In this section we will describe multi-particle simulation results. These will include the 
emittance growth, evolution of beam profiles, amplitude growth at different initial ampli- 
tudes, and also the growth of the variance in action at these initial amplitudes. This will 
allow us to probe both the macroscopic and microscopic beam behaviour. 

The simulations were performed with a simple numerical model consisting of six di- 
mensional linear transport between the two collision points, a sinusoidal longitudinal map 
through an rf cavity and weak-strong beam-beam interactions at the two IPs. The beam- 
beam interactions occur with a horizontal crossing angle at one IP and a vertical crossing 
angle at the second IP. The strong beam was assumed to have a Gaussian distribution in all 
three planes. Magnetic nonlinearities are not included, both to keep the model as simple as 
possible and also to avoid particle amplitudes from growing exponentially fast far from the 
beam core. Limiting amplitude growth to finite values allows us to keep all particles in the 
distribution and hence study the growth of the beam tails with good statistics. 
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Figure 1: (color) Emittance growth vs turns with tunes on the synchro-betatron resonances. 
Left: horizontal emittance growth, Right: vertical emittance growth. The power law fits 
and the exponents for the fits are also shown. 



3.1 Emittance growth and lifetimes 

Emittance growth was calculated by evolving ensembles of N particles (5000 <N< 20000) 
starting with Gaussian distributions in all planes. Typically 10,000 particles sufficed to 
obtain results that did not change much with a larger number of particles. The calculated 
emittance was the rms emittance, e.g. e x = [(x 2 )(x 2 ) — (xx 7 ) 2 ] 1 / 2 . Figure Q] shows the 
emittance growth with 20,000 particles on the two resonances. We find that the growth 
follows a simple power law, the fits are also shown in the figure. We observe that the 
horizontal emittance growth after 10 turns is more than 2.5 times larger on the 2(3 V x — 
2v s ) = 2 resonance than on the 4v x — 2v s = 1 resonance. The vertical emittance growth is 
much smaller than the horizontal, about a factor of five smaller for the first resonance and 
it is practically zero for the second resonance. 

By imposing a finite aperture restriction, we can find the escape time needed by particles 
to reach this aperture. This has been calculated for several different apertures and for both 
resonances. Apertures were placed from 5a to 10a at intervals of la. On the 2(3 v^ — 
2v s ) = 2 resonance, we find that about 7% of particles reach 8a, a handful reach 9a and 
none reach 10a. On the 4v x — 2v s = 1 resonance, about 4% of particles reach 6a, a few 
reach 7a and none reach 8a. The amplitude distribution of the particles reaching 8a on the 
first resonance and of the particles reaching 6 a on the second resonance are shown in the 
top plots of Fig[2j The initial distribution in each case was a Gaussian with 40,000 particles. 
On the 2(3 v x — 2v s ) = 2 resonance, the maximum of the amplitude distribution occurs close 
to 1.5a - an amplitude close to the lower edge of the resonance islands, shown later in Fig 
[TOl The minimum amplitude that reaches the aperture is 0.25 a. On the 4v x — 2v s = 1 
resonance, the corresponding peak in the amplitude distribution is close to 1.8a, also at the 
lower edge of the resonance islands seen in Fig. \TT\ The minimum amplitude that reaches 
the aperture on this resonance is 0.9a. 

The average escape time in the simulation may be interpreted as representing the beam 
lifetime. The bottom plot in Fig. [2] shows the average escape time (calculated with 40,000 
particles) as a function of the aperture amplitude for both resonances. The average es- 
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Figure 2: Top: Distribution of amplitudes reaching an aperture of 8c on the 2(3v x — 2v s ) = 
2 resonance (left) and an aperture of 6a on the 4v x — 2v s = 1 resonance (right). In both 
cases, the initial distribution was a Gaussian with 40,000 particles. Bottom (color): Average 
escape time for the two resonances at different apertures. 



cape time with 20,000 particles yielded similar values showing that these numbers have 
converged to stable values. The average escape time increases by an order of magni- 
tude or more for each increase in aperture by 1 a. The average escape time at 8 c on the 
2(3v x — 2v,y) = 2 resonance is about the same as at 6a on the 4v x — 2v s = 1 resonance. 
At a fixed aperture, the differences in escape times between the two resonances increases 
by about two orders of magnitude at 5 and 6 a and three orders of magnitude at 7a. One 
would expect this trend of increasing lifetimes to continue with higher order resonances. 



3.2 Beam profiles 



The beam profiles were found for the same distributions and resonances. The left plot 
in Fig [3] shows a mountain range view of the horizontal beam profiles (i.e. distribution 
function of the horizontal position), initially and then at other intervals up to 10 6 turns with 
tunes on the resonance 2(3v x — 2v$) = 2. After the initial time, the subsequent horizontal 
profiles develop long non-Gaussian tails which extend out to ±8a compared to the initial 
Gaussian distribution which was limited to ±3. 5a. The vertical beam profiles (not shown 
here) however stayed Gaussian and close to the initial distribution. The right plot in this 
figure shows the horizontal profiles but with tunes on resonance 4v x — 2v s = 1 . We observe 
that in this case as well that the tails are non-Gaussian and extend out to about ±6 a, not 
quite as far as on the first resonance. Again there is very little change in the vertical profile. 
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Figure 3: (color) Mountain range view of the horizontal beam profile initially and at subse- 
quent times. Particles were on the resonance 2(3v x — 2v s ) = 2 (left) and on the resonance 
4V X -2V S = 1 (right). 



We observe that the beam tails do not appear to change very much after the first 100,000 
turns or so. It is most likely that the regions of enhanced diffusion are depleted within 
these turns. The particles in the vicinity of the resonance islands are transported to larger 
amplitudes quickly and are detuned from the resonance. The amplitudes to which they 
move have much smaller diffusion, so the beam tails do not change much. As we will 
see in the next subsection, the evolution in the beam core shows growth even after several 
hundred thousand turns. However these particles do not migrate to the tails during the time 
duration followed. Thus we continue to observe emittance growth. 

In order to find distributions that can best fit the non-Gaussian tails, we first look to 
the Central Limit Theorem (CLT) which explains the ubiquity of the Gaussian distribution. 
This powerful theorem states that the distribution of a sum of a sequence of random, iden- 
tically distributed and independent variable with finite mean and second moment tends to 
a Gaussian distribution in the limit that the number in the sequence approaches infinity. 
Generalizing the CLT by dropping the requirement of a finite second moment leads to the 
family of Levy stable distributions [1201 . For applications in beam dynamics, these distri- 
butions will still have a finite second moment because they do not extend to infinity but are 
truncated at the beam pipe or the closest physical apertures. 

Levy stable distribution functions are defined by an inverse Fourier transform of a 
stretched exponentially decaying function in Fourier space 

1 f°° 

L a (z) = — exp[-izk-\k\ a ]dk, 0<a<2 (5) 

2% J-oo 

There is no known closed form expression for arbitrary values of a. Special cases include: 
the Lorentz distribution L\(z) while Li(z) is the Gaussian distribution. There are more 
general asymmetric versions of the Levy stable distribution with additional parameters but 
we shall not need them here. Some basic properties of these functions are fl2TJ 

• These functions are normalized : f " M dzL a (z) = 1 

• They are even functions : L a (— z) = L a (z) 
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Figure 4: (color) The final horizontal profile and a fit (blue) with a Levy stable distribution 



L a . Left: Resonance 2(3 V x 
a = 13. 



2y$) = 2 and a = 0.95. Right: Resonance 4v x - 2v s = 1 and 



At z = 0, L a (0) = i^r(^), which increases rapidly when a — > 0. 
At large values of z, the distributions decay as 



r 11 ,r(i + a) 

limL a (z) ~ — sin(-Tra) 

7T 2 



ll+a 



We find that the non-Gaussian horizontal profiles can be fit by these Levy stable dis- 
tributions L a . The left plot in Fig. |4] shows the fit of the final horizontal profile for the 
resonance 2(3v x — 2v s ) = 2 with a Levy stable distribution with parameter a = 0.95. This 
profile is narrower than a Lorentzian and decays at large x as |x| _L95 . The right plot in this 
figure shows the final distribution on the resonance 4v x — 2v s = 1 can also be fit by a Levy 
stable distribution with a larger central width and corresponding to a = 1.3. This profile is 
wider than a Lorentzian and decays at large x as |x| -2 3 . The Levy stable distributions were 
generated with a Mathematica package [fT9l . 

It is known G21 that the Levy stable distributions serve as Green's functions to frac- 
tional diffusion equations for a density p (jc, t) of the type 



(6) 



where % is a constant diffusion coefficient and ^D" is the Riemann-Liouville fractional 
space derivative of order a given by, 



l 



p(*0 



r(2-a) ^ 2 i-oc {x-x') a - 1 
The solution of the fractional diffusion equation above is 

/oo 
L a (z)p (x-(xt) l/a z)dz 
-co 



(7) 



(8) 



where po(x) is the initial density. Levy stable distributions have also been shown to be so- 
lutions of other fractional diffusion equations [23] . There is no reason to believe that either 
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Figure 5: (color) Plots of the initial (red) and final (blue) beam distributions in the hori- 
zontal plane at the resonance 2(3v x — 2v s ) = 2 (left) and resonance 4v x — 2v s = 1 (right) 
respectively. Initially 4000 particles each were placed at horizontal amplitudes from 0.5 
to 4a in steps of 0.5a. The vertical amplitude was kept constant at 0.1a. In the plots the 
vertical scale has been truncated to 400 in order to show clearly the particle numbers in the 
final distribution at large amplitudes. 



Equation © or of the type in reference [|23l are appropriate for our problem. However the 
fact that the long time beam profiles are described by these Levy distributions is our first 
indication that the amplitude growth process may be described by an appropriate fractional 
diffusion equation rather than the regular diffusion equation. In Appendix A we derive a 
different fractional diffusion equation that may describe the dynamics observed here. 



3.3 Growth at individual amplitudes 

We now take a closer look inside the beam distribution to determine how the amplitude 
growth changes with amplitude. Instead of a Gaussian distribution in phase space, we 
consider delta function distributions in action. We select a discrete number of horizontal 
actions and at each action we place 4000 particles uniformly distributed in angle. The 
vertical amplitude was kept constant at 0.1a for all particles. The initial distribution in 
transverse action angle space can be written as 

P (/„ 6 X , Jy, By) = 8(jy- J O .l)P(e x )P(0y) £ 5 (J x ~ /,) (9) 

i 

where Jq_\ is the action at an amplitude of 0.1a, P(6 X ) is a uniform distribution in the 
horizontal angles etc. The initial longitudinal variables were chosen to be the same for all 
particles: z = la s , 8p/p = lo p in these simulations. We let these distributions evolve and 
record the final distribution in amplitude after 10 6 turns. The left plot in Fig |5] shows the 
initial (red) and final (blue) distributions for resonance I. We observe that particles at 0.5a 
stay close to their initial amplitude. At la, many particles have moved to larger amplitudes 
but a sizable fraction stay in their original neighbourhood. This shows a large variation 
in final amplitude depending on their initial angle or sensitivity to their initial conditions. 
It suggests that motion in the neighbourhood of la could correspond to bounded chaos. 
At amplitudes of 1.5a and higher, the vast majority of particles have migrated to larger 
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Figure 6: (color) Variance in the horizontal actions over time at tunes corresponding to the 
resonances 2(3v x — 2v s ) = 2(left) and 4v x — 2v s = 1 (right). Also shown are monomial fits 
to the data. Note that the variances are plotted on a logarithmic scale. All variances are 
zero initially but the zero is suppressed here. 

amplitudes up to 8c and depleted the initially populated regions. There is a broad local 
maxima in the final distributions at ~ 7a. The right plot in Fig. [5]shows the corresponding 
results for resonance II. The results are qualitatively similar with some differences. The 
initial amplitude with large variation in final amplitude is closer to 2a and the largest 
amplitude reached is about 7a. On this resonance there remain local spikes at 2.5 and 3a 
showing that diffusion at these amplitudes is weaker than in the first resonance. 

3.4 Variance of the action and diffusion type 

We now examine the diffusion from individual amplitudes. In regular diffusion the variance 
of the diffusing quantity, here the action, grows linearly with time which allows one to 
define time independent diffusion coefficients D(J) = ((AJ) 2 )/At. We check the validity 
of this assumption for the beam-beam driven SBRs. using the same initial distributions as 
used in Fig [5J Variances are calculated over particles at the same initial action. Figure [6] 
shows the growth in the variance of the horizontal action at several initial actions for both 
resonances. The vertical amplitude was constant at y = 0.1a. Initially the variance is zero 
at all actions but then grows at different rates depending on the action. The growth in the 
variance is not linear at any action. In most cases there is a sharp initial transient growth 
which is followed by a slower long term growth. This long term growth can be modeled 
(again in most cases) by a power law behavior of the form 

((AJ x ) 2 }~C x t P % ({&J y ) 2 )~Cyt P y (10) 

where the coefficients (C x ,C y ) and the powers (p x ,P y ) depend on the initial action. Ex- 
ponents less than 1 indicate sub-diffusive behavior while exponents greater than 1 imply 
super-diffusive motion Figure |6] also shows the fits with this power law. Growth of the 
variance in the vertical action can also be fit by a single power law with small values of 
(Cy,p y ) showing that there is no appreciable diffusion in that plane. On the resonance 
2(3v x — 2v s ) = 2, there is significant growth in the action at amplitudes of 2 and 2.5a com- 
pared to neighboring actions both lower and higher. The exception to the single power law 
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fit occurs at x = 1 a where the variance stays nearly constant after the initial transient and 
then after about 400,000 turns grows by an order of magnitude over the next 600,000 turns 
but with oscillations in the variance. These oscillations occur because of the large sensi- 
tivity to the initial angle at this amplitude. The oscillations decrease significantly when 
the number of particles at the same initial action is increased from 4000 to 20000 parti- 
cles, which results in a more complete sampling of the initial angle. Simulations show 
that this greater sensitivity to the initial angle is also present at amplitudes in the range 
1.0a < x < 1.3a with y = 0.1a. The fits to a power law in this zone are applied after the 
variance starts to grow rapidly but with 20000 particles. The average action with initial 
\x\ = 1.0a grows about 10% after 10 6 turns while the average action with initial \x\ = 1.5a 
grows by about a factor of two over this time. So the narrow zone around \x\ = 1.0a 
corresponds to a zone of bounded chaos. 

At the resonance 4v x — 2v s = 1, the growth in variance is largest in the range x = 
2.5 — 3a and drops for both smaller and larger initial actions. The large oscillations in the 
variance occur in a range around x = 2.0a and again these oscillations are reduced when 
the number of particles is increased from 4000 to 20,000. For this resonance, the zone 
around \x\ = 2.0a is a zone of bounded chaos. Similar behaviour is seen at other values of 
y but the width of the zone of bounded chaos changes. 

The exponents in the power laws were calculated for several values of the horizontal 
amplitude and for different vertical initial amplitudes. Fit [7] shows the exponents for both 
resonances. On the resonance 2(3v x — 2v s ) = 2 there is a spike in the exponent to values 
well above 1 in the regions of bounded chaos for y = 0.1,0.5a suggesting super-diffusive 
behavior. Above the zone of bounded chaos, the exponent falls well below 1 suggesting 
sub-diffusive behavior. At y = la the exponent stays well below 1 for all x showing that 
zones of bounded chaos have disappeared. On the 4v x — 2v s = 1 resonance, the exponent 
rises above lonly in a narrow zone around x = 2a at y = 0.1a. At y = 0.5a the exponents 
stay well below 1 at all x with a small spike at x = 2a. The motion is sub-diffusive at all 
x values studied when y = la. Since the super-diffusive regions are narrow, it is possible 
that they may appear for \y \ > 1 a when the motion is studied with a finer resolution or even 
when the longitudinal variables are changed. We remark that we have observed here three 
different signatures of bounded chaos: large variations in final amplitude when starting 
from the same initial amplitude (seen in Fig [5]), large oscillations in the action variance 
over time (seen in Fig. |6]and a spike in the power law for the growth of the variance (seen 
in Fig. |7]). These signatures apply to an ensemble of particles at the same amplitude but 
different initial angle as opposed to the Lyapunov exponent criterion which is applied to a 
pair of particles that are initially infinitesimally close. 

The picture that emerges is that near synchro-betatron resonances, phase space is di- 
vided into several zones. At small amplitudes there is no diffusion. At larger amplitudes 
there is a zone of bounded chaos with super-diffusive motion. The next zone outward in 
phase space is wider with sub-diffusive motion. Finally at even larger amplitudes, the mo- 
tion becomes linear again and consequently there is no diffusion. Fig. [8] shows a qualitative 
sketch of these different zones. The width of the super-diffusive zone with bounded chaos 
depends on the resonance, on the amplitude of the orthogonal transverse amplitude and on 
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Figure 7: (color) Exponent p x of time in the power law fits of the horizontal action variance 
vs the initial horizontal amplitude for different initial values of the vertical amplitude. The 
left figure corresponds to the resonances 2(3v x — 2v s ) = 2 and the right to the resonance 
4v x — 2v s = 1. The exponent spikes above 1 in a very narrow range of horizontal ampli- 
tudes. Exponent values below 1 indicate sub-diffusive behavior while those above indicate 
super-diffusive behaviour. 

the values of the longitudinal variables. 

The fact that the sub-diffusive regions seem to be dominant in this perturbed Hamilto- 
nian system should not be unexpected due to the existence of hyperbolic fixed points and 
the existence of perturbed KAM tori. These fixed points and tori lead to orbits which stay in 
their vicinity for long time periods and consequently to slower growth. Similar phenomena 
have been reported for the standard map by Balescu [24] . 

4 Statistics of single particle behavior 

We saw in the previous section that in most regions of phase space, the variance grows 
slower than linearly with time. If we define an instantaneous or 'running' diffusion coef- 
ficient [|24l as Dj x = (1/2) (9 (AJ%) /dt, then this coefficient would be time dependent and 
would vanish in the very long time limit. Near both resonances we did not observe any 
zone of regular diffusion with constant diffusion coefficients. We also saw that the beam 
profile was given by a Levy stable distribution which is known to be the solution of a frac- 
tional diffusion equation. These suggest that the dynamics near these resonances cannot 
be described by the regular diffusion equation but instead that the diffusion is anomalous 
which needs a different diffusion equation. In order to test this possibility in more detail, 
we will examine the validity of the assumptions behind the regular diffusion equation. 

4.1 Continuous Time Random Walks 

The regular diffusion equation arises after assuming that the particle dynamics can be mod- 
eled as a classical random walk following a Markov process. This implies that particle 
jumps occur at regular time intervals and there is a well defined time scale such that events 
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Figure 8: (color) Qualitative sketch of phase space divided into zones of no diffusion, super- 
diffusion and sub-diffusion when the dynamics is dominated by a beam-beam synchro- 
betatron resonance. 

separated in time by longer than this time scale are uncorrelated. It then follows that the 
particle density is governed by the well known Chapman-Kolmogorov master equation. 
From this master equation and a few more assumptions (e.g. on the smallness of the dis- 
placements etc.) the regular diffusion equation follows. See Appendix A for a sketch of 
this derivation. 

A well known alternative to the standard random walk picture is the Continuous Time 
Random Walk (CTRW) model introduced by Montroll and Weiss [25J to consider processes 
where both the times at which jumps occur as well as the sizes of the jumps in space are 
random functions. A review of CTRW and connections to fractional diffusion equations 
may be found in ll26l . 

A general dynamical process may not have a characteristic time scale. In those cases a 
Markov description may not be applicable. The CTRW model introduces the concepts of a 
probability distribution w for the waiting times before a jump occurs and a probability dis- 
tribution *F for the size of a jump In beam dynamics there is no diffusion when the motion 
is linear and the usual Courant-Snyder actions are conserved. Consequently it makes sense 
to define the jumps in action space when the motion is nonlinear and diffusive. Hence we 
define w(t, J) At to be the probability that a particle waits for a time between t and t + At at 
action J before making a jump, and define ^(AJ; J,f)AJ to be the probability of making a 
jump by AJ at the action J at time t. These distributions are normalized, i.e. 



The concept of a waiting time endows the system with memory. The CTRW model reduces 
to the classical random walk model on which the regular diffusion equation is based, when 



These waiting time and jump size distributions can be used in many cases to determine 




(11) 



the waiting time follows an exponential behavior in time e 
scale T. 



///t with a characteristic time 
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Turns [xlOOO] Turns [xlOOO] Turns [xlOOO] 



Figure 9: Time series of the horizontal amplitude with x = 1.5a, y = 1.5c Left plot: We 
see the first large jump in the amplitude from ~ 1.5c to ~ 6.5a occurs after 140,000 turns. 
Middle plot: This zoomed in plot shows sequences of jumps from small to large amplitudes 
and back interspersed with intermittent periods of small amplitude changes. Right: The last 
200,000 turns during an evolution over 2 x 10 6 turns. Here we see excursions between 1 to 
8 a. 

the evolution followed by the density distribution p(J,f). The canonical CTRW model 
assumes a power law waiting time distribution, a Gaussian for the jump size distribution 
and a constant diffusion coefficient. These lead to a fractional diffusion equation for the 
density ll2~6ll . In our case the dynamics near resonances is sufficiently complicated that we 
need to establish the evolution equation for the density from first principles. We therefore 
need to determine the forms of the jump size distribution and the waiting time distributions 
from the dynamics. Simulations discussed in the rest of this section are used to extract 
these distributions. 

A check of the whether the CTRW model may be applicable here can be done by ex- 
amining the time series of single particles. Fig |9] shows one example of a time series of the 
amplitude a/2/3 x 7 x for a single particle on the resonance 2(3 V x = 2v s ) = 2. The left plot 
shows that a particle may perform small amplitude quasi-periodic oscillations for a while 
before a major qualitative change occurs. The middle and the right plots show that step 
sizes can be large (several a), of varying amplitude, and there are intermittent sequences 
of varying duration where there are smaller steps. The time dependent behaviour of this 
sequence and the non-locality of the changes establish that this is a process with a distri- 
bution of waiting times and a distribution of action step sizes, the key ingredients of the 
CTRW model. 

4.2 Jump size distributions 

We now calculate the jump size distributions by following a single particle for 10 6 turns 
and find the changes Ax, AJ X in position and action per turn. Fig. [T0l shows the phase space, 
and jump distributions of Ax, AJ X on the resonance 2(3v x — 2v,y) = 2 with initial values of 
x = (0.2,2, 8)a. At the smallest initial position xq = 0.2a, the phase space is a distorted 
ellipse with no trace of the resonance island; motion here is quasi-linear. The plot for the 
distribution function of Ax also has the distribution function for a periodic function shown 
in dotted lines. When the argument of a periodic function like sine or cosine is sampled 
from a random distribution, the distribution function for the periodic function / has the 
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Figure 10: Phase space (left), jump size distribution in x (middle) and jump size distribution 
in action J x (right) on the resonance 2(3v x — 2v s ) = 2. The initial value of x changes going 
from top to bottom as x = (0.2,2.0, 8.0) a. The initial value of y = 0.1a is the same in all 
these plots. The distribution in AJ X is plotted on a semi-log scale and the abscissa is in units 
of AJ x /J<j where J a is the action at 1 o x . 



form 

p(f) ~ -7=1=, \f\ < 1 (12) 
V 1 -/ 

The distribution function has local maxima wherever the function itself becomes stationary, 
so that many more points are sampled from the neighbourhood of these stationary points. 
Since the motion at small and large amplitudes is quasi-periodic in our model, it is to be 
expected that the distribution in Ax is close to that of a periodic function. The distribution 
function for AJ X is plotted on a semi-log scale and shown as discrete points, for greater 
clarity. At xq = 0.2a, the distribution for AJ X lies on a single curve but not given by any 
simple expression. As the particle's initial position increases to 2 a, the nonlinearity of 
the beam-beam force manifests and we see resonance islands in phase space and large 
excursions. The distribution function for Ax undergoes a qualitative change to resembling 
a parabolic curve but with a dip in the center and with peaks close to the center. The 
distribution function for AJ X now falls on two separate curves. Similar distributions for 
Ax, AJ X are seen for initial particle amplitudes in the range 1 .5a < |jco | < 6.5a. At xq = 8a, 
the phase space returns to a distorted ellipse with considerable smear, and the distribution 
functions also resemble those seen at xq = 0.2a. 

Fig. [EH shows similar plots on the resonance 4v x — 2v s = 1 with initial values of x = 
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Figure 11: Similar plots as in Fig [TO] but on the resonance 4v x — 2v s = 1. From top to 
bottom, the different initial values of x = (0.5, 5, 8)<7. The initial value of y = 0. Id is the 
same in all cases. 



(0.5,3.5, 8)<7. Again, we see a qualitative change in the distribution functions when the 
motion is strongly nonlinear in the presence of the resonance islands. The shapes of the 
distributions in Ax are similar to those seen for the previous resonance and the distribution 
of AJ X also lies on separate curves at intermediate amplitudes. These suggest that there is a 
universal character to the jump distributions which mirrors the behavior in phase space. 



4.3 Waiting time distributions 

The waiting time distribution is the important distribution that determines the nature of the 
diffusion process. As remarked earlier, a waiting time distribution that follows an exponen- 
tial law reduces to a Markov process, otherwise the process is non-Markovian. The waiting 
time for each initial amplitude is found here by tracking a particle at that amplitude for 10 6 
turns. The phase space region in action angle coordinates that is visited by the particle is 
divided into different zones and the time that the particle stays in the zone before leaving 
is one instance of the waiting time. The choice of the width of the zone is somewhat arbi- 
trary since there is no dynamics dependent action scale which is applicable to all of phase 
space. For example, the resonance width is not relevant at small or large amplitudes and if 
there were multiple resonances, there would be multiple widths. We therefore calculated 
the waiting time distribution twice, once with a chosen width such that there was enough 
statistics in each zone and the second time with twice the width. In most cases we found 
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Waiting time [turns] Waiting time [turns] 

Figure 12: (color) Waiting time distribution (on a log-log scale) for the 2(3v x — 2v s ) = 2 
resonance (left) and 4v x — 2v s = 1 resonance (right). The distributions are calculated at 
amplitudes where there is significant diffusion of particles to larger amplitudes. 



that the parameters of the distribution change by less than 10%; we take this to be a sign of 
convergence of the distribution. We find that the exponential function is not a good fit to the 
distribution for either resonance. The results for a fit to a power law distribution are shown 
in Figure \Y2\ The distributions are plotted on a log-log scale for several initial amplitudes 
where there is significant amplitude growth. On the 2(3v x — 2v s ) = 2 resonance, most of 
the points (with the exception of the single occurrence events with long waiting times) lie 
on straight lines showing that a power law is a reasonable fit. The power law exponents for 
the different amplitudes are close. For the amplitudes shown in this figure, the waiting law 
distributions are 

w(0~r a , 2.4<a<2.7 (13) 

On the 4v x — 2v s = 1 resonance, the waiting time distribution can also be fit by a power 
law distribution but the range of variation in the exponent a is larger: 1.4 < a < 2.7. The 
greater variability in the exponent is expected to have an impact of the dependence of the 
diffusion rate at different amplitudes on this resonance. 



5 Fractional diffusion equation 

Since the waiting time distribution suggests that the transport near resonances is non- 
Markovian, we need to establish an alternative to the regular diffusion equation. For a 
Markov process, the regular diffusion equation is obtained from the Chapman-Kolmogorov 
master equation, a derivation is sketched in Appendix A. In Appendix A we also derive 
a different master equation using general jump size and waiting time distributions for a 
CTRW process in action space following a method outlined in [|27l . The master equation 
for the density in action angle space that we obtain is 

d 1 f f 1 

p(j 5 0)=_ / / dAJdAe x ¥(J-AJ,0-Ae;AJ,A0)L t p(J-AJ,0-Ae,t)--L t p(J,e,t) 
at x J J T 
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where L t is an integral operator given by 




sw(s; J, 9) 



P(J,M 



(15) 



l-w(s;J,6) 



Here _1 is an inverse Laplace transform, T is a time parameter in the waiting time dis- 
tribution w(t, J), h^s; J),p(s; J) are the Laplace transforms of the waiting time distribution 
and the density respectively. In the Appendix we then show that expanding this master 
equation in a Taylor series in the same manner as is done for the Chapman-Kolmogorov 
equation, the following fractional diffusion equation is obtained for a power law waiting 
time distribution w(t, J) ~ t~ a ^> 



Here the exponent a depends on the action J which will be true in general and D^i are 
action dependent diffusion coefficients, defined in the appendix. It remains to be verified 
that this fractional diffusion equation describes the dynamics near resonances, as seen in 
the particle tracking simulations. However, this diffusion equation has been derived under 
general considerations of a CTRW process which the dynamics near the SBR resonance 
appears to follow. Given the large variations in the diffusion coefficients, the solution 
of this diffusion equation will likely require a special purpose numerical algorithm. The 
density can then be used to calculate the beam lifetime and various moments such as the 
emittance . 

6 Discussion 

We have studied the detailed transport process near two low order horizontal synchro- 
betatron resonances driven by beam-beam interactions at a crossing angle. We found that 
the horizontal beam profiles develop long beam tails. The horizontal beam distribution 
evolves from an initially Gaussian distribution to a Levy stable distributions on both reso- 
nances. The Levy stable distributions are solutions of simple fractional diffusion equations 
which describe some anomalous diffusion processes. The evolution of the variance in ac- 
tion at several initial values characterizes the nature of the diffusion in phase space. At 
small amplitudes there is no diffusion, then there is a narrow region where the motion is 
super-diffusive (the variance grows faster than linearly with time), followed by a broad 
region where the motion is sub-diffusive (the variance grows slower than linearly with 
time) and finally no diffusion at large amplitudes. The width and the location of the super- 
diffusive region depends on the resonance, the width is narrower for the weaker resonance. 
This super-diffusive region is also marked by signatures of bounded chaos and particles 
do not experience large amplitude growth. For both resonances, this region is located at 
the lower edge of the resonance islands. The broad sub-diffusive region abuts the super- 
diffusive region and continues until about 5-6 a depending on the resonance. Here particles 
do migrate to larger amplitudes. We do not observe regular diffusion anywhere in phase 
space on either resonance with the particle distributions we used. 




(16) 
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The jump size distribution and the waiting time distribution, key ingredients of a con- 
tinuous time random walk process, were found by analysis of single particle tracking data. 
The jump size distributions for both resonances were similar - in the linear regions of phase 
space, the distributions in Ax are close to the arcsine distribution while in the nonlinear re- 
gions they have a more complex shape. The similarity of these distributions for the two 
resonances suggests that these may be universal features near such resonances. When the 
waiting time distributions follows an exponential law, the stochastic process is Markovian. 
We find that the waiting time distribution follows instead a power law, again for both res- 
onances. Since the process is non-Markovian, the regular diffusion equation cannot be 
used to describe the evolution of the density. For a general CTRW process, we derived a 
master equation in action-angle space which is applicable to processes with arbitrary jump 
size and waiting time distributions. A fractional diffusion equation was derived from this 
master equation. Numerical solutions of this diffusion equation will allow computations of 
beam observables such as lifetimes and emittance growth. 

This model can be tested against beam observations when anomalous diffusion is sus- 
pected. Comparison of beam profiles with Levy stable distributions would be a first check. 
Another indicator would be if the emittance of pencil beams grow nonlinearly with time. 
This could then be followed by measurements of diffusion coefficients at different am- 
plitudes, using them in the fractional diffusion equation and comparing the numerically 
calculated emittance growth and beam lifetime with the measured values. 

In this article we considered low order synchro-betatron resonances so as to observe 
effects on a short time scale. Based on comparisons of the two resonances studied here, we 
expect that the physics at high order resonances (and hence more applicable to operational 
accelerators) will be similar but on longer time scales. When multiple such resonances 
are present simultaneously, the diffusion is likely to be anomalous but the phase space 
dynamics will be more complicated. It is possible that the physics near space charge driven 
resonances may be similar to that obtained here but that remains to be investigated. 
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A Appendix: Regular and fractional diffusion equations 



We briefly summarize the derivation of the diffusion equation in action-angle space. We 
assume a Hamiltonian description H(J,9) which has been perturbed from an integrable 
Hamiltonian Hq[S). Let ^(J, 0; AJ, AO be the transition probability for the action-angle 
variables to change from (J, 0) to (J + AJ, + AO) in time At. The first major assumption 
is that the dynamics is Markovian. For a Markov process, the particle density distribution 
at time t + At only depends on its instantaneous state at t and is independent of its previous 
history provided At is longer than a characteristic time 1. Under this assumption, the den- 
sity p ( J, , t) at time t + At can be found by summing over all possible transitions in time 
At. This results in the Chapman-Kolmogorov equation for the density 

p(J,0,t + At) = j J p(J-AJ,0-A0,t)*¥(J-AJ,0-A0;AJ,A0)d(AJ)d(A0) (A.l) 

Here *F is the transition probability of jumps (AJ,A0). Further assumptions need to be 
made including i)the angles evolve on a faster time scale than the actions and their corre- 
lation decays rapidly, ii) the density in the long time limit is independent of the angle iii) 
the transition probability can be factorized in the form *P( J, ; AJ, AO) = x ¥j(J; AJ) S (AO — 
OAt) iv) the changes in action and angle AJ,A0 are small during a time interval At. Ex- 
panding the LHS and the RHS of Equation (IA. II) . keeping up to second order terms and 
then taking the limit At — > 0, we obtain the Fokker-Planck equation 

| = -V,.[Ap] + ££^1^] (A.2) 

where the drift A and diffusion coefficients D are defined as 

A ^= S n %T> D M = AT lim (AJ>= /AJ^ y (J;AJ)JJ 

At AJ->-0,A?-s-0 2 At J 

(A.3) 

Here At is understood as a time shorter than a time scale over which the density distribution 
evolves but longer than the time over which angle correlations decay. 

For Hamiltonian systems, there is a relation between the drift coefficient and the diffu- 
sion coefficients If28llj9ll 

then the Fokker-Planck equation simplifies to the diffusion equation 

The assumptions of Markovian behavior and the smallness of the changes in action-angle 
variables are crucial for the validity of this regular diffusion equation. If these assumptions 
are invalid, then this diffusion equation may not be the right model for the density evolution. 
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We now consider a more general master equation for a CTRW process in action angle 
space with arbitrary jump size and waiting time distributions. We use a method outlined 
in [|27l . It uses two basic balance conditions: the first states that a change of density arises 
from the difference in the incoming flux r + (J, 9,t) and the outgoing flux r~(J, Q,t). 

§- f p(3,o) = r+{j,e,t)-r-{j,e,t) (A.6) 

The second balance condition states that the influx is composed of the outflux of particles 
from all other phase space locations to that location 

r+(j,e,*) = y ydAjdA0¥(j-Aj,e-Ae;Aj,A0)r-(j-Aj,e-Ae,o (A.7) 

The outflux at (J, 9,t) has contributions from particles that were present initially but left 
after waiting for time t and those that arrived later before leaving 



r-(j,e,0 



w(f;j,e)p(j,e,o)+ [ \{t-t' 'jj.ejr 4 (j,e, t')dt' (A.8) 

Jo 



Substituting Eq. (IA.8I) in Eq. (|A.6I) and taking the Laplace transform, we obtain for the 
outflux 

1 



r-(j,e,t) = & 



-i 



sw(s-,j,e) 

-p{j,e,s) 



= -Ltp(J,d;t) (A.9) 

X 



l-w(s;J,6) 

Here w(s;J, 9) and p(J, 9,s) are the Laplace transforms in s space, x is a relevant time 
parameter in the waiting time distribution, and Jzf _1 is the inverse Laplace transform. The 
last equality in this equation defines the integral operator L t . Substituting this back in Eq. 
(IA.6I) and using Eq. (IA.7l) we obtain 



d_ 

dt 



p(J,0)=Iy J JAJJA0^(J-AJ,0-A0;AJ,A0)L f p(J-AJ,0-A0,O-^p(J,0,O 

(A.10) 

This is the modified master equation for the density. 



Now we derive the modified diffusion equation from this master equation. We expand 
the RHS of Eq. dA.10l) in a Taylor series and keep up to second order terms. As before we 
define the coefficients 

A(J)=BmM, %(J)= h !« (A. 1 1) 

AJ^O, X AJ^O 2 X 

We assume that the same relation as in Eq. (|A.4b between the drift and diffusion coefficients 
holds. Then we have as the modified diffusion equation 

In cases where (1/ x)Ltp = p, this is the regular diffusion equation. 
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Consider now two examples of a waiting time distribution, first an exponential waiting 
time 

w(r) = I=p[-i], ^*( S ) = I(- ? L_) (A.13) 
The integral operator simplifies to 

-Up = 5T 1 [y^yP (J, s)} = p ( J, t) (A. 14) 

i.e. the modified diffusion equation reduces to the regular diffusion equation. 
Now consider a power law waiting time 

W (fj) = -(-r a{3) (A.i5) 

T T 

Here we let the exponent a be action dependent. In the long time limit t — > °° or equiva- 
lently s ->■ 0, 

-L f p=r(l-a(J))^- 1 [^ (J) P(J^)]=oA a(J) p(J,0 (A.16) 

T 

Here T is the Gamma function and qD"^ is a Riemann-Liouville fractional derivative in 
time defined below. The diffusion equation for p is 



t=IE^[-4]o^P(J,o 



* ' * * 1 r " ^//' piJ : //, .l (A.i7) 



r rr dJ k K wr(i-a(j)) 



dtJo (t-t') a{3) 



This is a non-local in time (due to the waiting time distribution) integro-differential diffu- 
sion equation for the density. 
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